/** This program estimate a factor model using a set of measurement equations. **/


void getpar_loglike(const gsl_vector* x, const int n){
	
    int ivec = 0;	
	
    	
	int NUM_measures_control_now = NUM_measures_control; 
	int im_measure_control_now   = 0; 
	
	int icount =0; 
	
	
	if (NUM_icat_unobs != 1) {printf("\n\n Errors in getpar_loglike(): NUM_icat_unobs = %d != 5 \n\n exiting \n ", NUM_icat_unobs); exit(1); }

	
	double varunobs[NUM_unobs_loglike]; 
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++) {
		varunobs[iunobs] = exp( gsl_vector_get( x, ivec  ) ); 			ivec++; 
    }
	
	if (varunobs[1]>1.0) varunobs[1]=1.0;

	double temp1 = ( gsl_vector_get( x, ivec  ) ); 			ivec++;
	clcov_unobs[0] = varunobs[1] * ( 2.0 *temp1 /(1.0+temp1 ) - 1.0); 
	
	{
		int iunobs; 

		iunobs = 0; 
		cmu_unobs[0][iunobs] =  exp( gsl_vector_get( x, ivec  ) ); 		ivec++;   
		
		
		iunobs = 1; 
		cmu_unobs[0][iunobs] = exp( gsl_vector_get( x, ivec  ) ); 		ivec++;   

	}
	


	calpha_measures_thetac[0] = 1.0; 
	
	for(int im = 1; im < NUM_measures_thetac; im++){
	
		calpha_measures_thetac[im] = gsl_vector_get( x, ivec  ); 		ivec++;   
	}
	

	for(int im = 0; im < NUM_measures_thetac; im++){
		
		cmu_measures_thetac[im] = gsl_vector_get( x, ivec  ); 		ivec++;   
		
		csigma_error_thetac[im] = exp( gsl_vector_get( x, ivec  )); 		ivec++; 	// measurement errors

		
		// covariance of thetac 
		for (int ix =im_measure_control_now; ix < NUM_measures_control_now; ix++){
			cbeta_measures_thetac[im][ix] = gsl_vector_get( x, ivec  ); 	ivec++;   
		}
		
		
	}



	calpha_measures_thetan[0] = -1.0; 

	
	for(int im = 1; im < NUM_measures_thetan; im++){

		calpha_measures_thetan[im] = -exp( gsl_vector_get( x, ivec  ) ); 		ivec++;   
		
		if (calpha_measures_thetan[im] < -2.0) calpha_measures_thetan[im] = -2.0; 
		
	}

	
	for(int im = 0; im < 2; im++){
		
		csigma_error_thetan[im] = 1.0; 
	}
	for(int im = 2; im < NUM_measures_thetan; im++){
		
		csigma_error_thetan[im] = 1.0; 
		
		cutoff_measure_thetan[im-2] = exp ( -gsl_vector_get( x, ivec ) );     ivec++;
		if (cutoff_measure_thetan[im-2] > 3.0) cutoff_measure_thetan[im-2] = 3.0; 
	}
	
	
	for(int im = 0; im < NUM_measures_thetan; im++){

		cmu_measures_thetan[im] = gsl_vector_get( x, ivec  ); 			ivec++;   
			
		// covariance of thetan 
		for (int ix =im_measure_control_now; ix < NUM_measures_control_now; ix++){
			cbeta_measures_thetan[im][ix] = ( gsl_vector_get( x, ivec ) ); 	ivec++;   
		}

	}
	
	
		
#ifdef DEBUG
    if (ivec != n){
	std::cerr << "\n\n Error in getpar_loglike(): Nvar_loglike = " << n << " ivec = " << ivec <<" \n exiting ...\n"; 
	exit(1); 
    }
#endif 
	
	
	
	for (int iunobs =0 ; iunobs < NUM_unobs_loglike; iunobs++){
		cmu_unobs0[iunobs] = 0.0; 			
		for (int icat=0; icat < NUM_icat_unobs; icat++) cmu_unobs0[iunobs] -= cmu_unobs[icat][iunobs]*DATA_init_icat_control_mean[icat]; 
	}
		
	
	clvar_unobs[0] = sqrt( varunobs[0] ); 	
	clvar_unobs[1] = sqrt( varunobs[1] - clcov_unobs[0]*clcov_unobs[0] ); 	
			
	//... mat_L is the lower triangular matrix with real and positive diagonal entries from Cholesky decomposition of Variance matrix: Sigma = L * L'
	gsl_matrix_set_zero (MAT_lcov_unobs); 	
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
		gsl_matrix_set( MAT_lcov_unobs, iunobs, iunobs, clvar_unobs[iunobs]); 
		
	}
	icount =0; 
	for (int iunobs=1; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol =0; icol < iunobs ; icol++){
			gsl_matrix_set( MAT_lcov_unobs, iunobs, icol, clcov_unobs[icount]); 
			icount++; 
		}
	}
	
	gsl_matrix_set_zero (MAT_cov_unobs);
	gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, MAT_lcov_unobs, MAT_lcov_unobs, 0.0, MAT_cov_unobs);

	
}


void getvec_loglike(gsl_vector* x, const int n){

    int ivec = 0;	
    		
	int NUM_measures_control_now = NUM_measures_control; 
	int im_measure_control_now   = 0; 

	


	double varunobs[NUM_unobs_loglike]; 
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++) {
        varunobs[iunobs] = gsl_matrix_get(MAT_cov_unobs, iunobs, iunobs);
        gsl_vector_set( x, ivec,  log( varunobs[iunobs] ) );  		ivec++;
    }
    

	double y1 = clcov_unobs[0]/varunobs[1] + 1.0; 
	gsl_vector_set( x, ivec,  ( y1/(2.0-y1) ) );  		ivec++;


	{
	
	int iunobs; 
	
	iunobs = 0; 

	gsl_vector_set( x, ivec, log( cmu_unobs[0][iunobs]) ); 		ivec++;   
		
	iunobs = 1; 
	gsl_vector_set( x, ivec, log( cmu_unobs[0][iunobs]) ); 		ivec++;   

	}


	
	for(int im = 1; im < NUM_measures_thetac; im++){

		gsl_vector_set( x, ivec,  calpha_measures_thetac[im]);  		ivec++;   
	}
		

	
	for(int im = 0; im < NUM_measures_thetac; im++){
		
		gsl_vector_set( x, ivec,  cmu_measures_thetac[im]);  		ivec++;   
		
		gsl_vector_set( x, ivec,  log(csigma_error_thetac[im]));  		ivec++; 	// measurement errors

		// covariance of thetac 
		for (int ix =im_measure_control_now; ix < NUM_measures_control_now; ix++){
			gsl_vector_set( x, ivec,  cbeta_measures_thetac[im][ix]); 	ivec++;   
		}
			
	}




	 
	for(int im = 1; im < NUM_measures_thetan; im++){

		gsl_vector_set( x, ivec,  log(-calpha_measures_thetan[im] ));  		ivec++;   
	}


	for(int im = 2; im < NUM_measures_thetan; im++){
		
		gsl_vector_set( x, ivec,  -log(cutoff_measure_thetan[im-2]));  		ivec++; 	
	
	}
	
	
	for(int im = 0; im < NUM_measures_thetan; im++){

		gsl_vector_set( x, ivec,  cmu_measures_thetan[im] );  			ivec++;   
		
	
		// covariance of thetan 
		for (int ix =im_measure_control_now; ix < NUM_measures_control_now; ix++){
			gsl_vector_set( x, ivec,  (cbeta_measures_thetan[im][ix]) ); 	ivec++;   
		}
	
	}

	
#ifdef DEBUG
    if (ivec != n){
	std::cerr << "\n\n Error in getvec_loglike(): Nvar_loglike = " << n << " ivec = " << ivec <<" \n exiting ...\n"; 
	exit(1); 
    }
#endif 



}


void readpar_loglike(const int n, ifstream &infile){

	double temp_x[n]; 

    for(int i=0; i< n && !infile.eof(); i++) infile >> temp_x[i];
        
    int i = 0;
 
	for(int im = 2; im < NUM_measures_thetan; im++){
 	cutoff_measure_thetan[im-2] = temp_x[i];		i++; 
 	}
 	
 	// unconditional distribution of unobservables
		gsl_matrix_set_zero (MAT_cov_unobs); 	
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
			gsl_matrix_set( MAT_cov_unobs, iunobs, iunobs, temp_x[i]); 		
			i++; 
		}
		int icount =0; 
		for (int iunobs=1; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol =0; icol < iunobs ; icol++){
			gsl_matrix_set( MAT_cov_unobs, iunobs, icol, temp_x[i]); 
			gsl_matrix_set( MAT_cov_unobs, icol, iunobs, temp_x[i]); 
			icount++; 
			i++; 

		}}	
		
		gsl_matrix_memcpy (MAT_lcov_unobs, MAT_cov_unobs); 
		gsl_linalg_cholesky_decomp (MAT_lcov_unobs); 
	for (int iunobs=0; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol = NUM_unobs_loglike-1; icol > iunobs ; icol--){
			gsl_matrix_set( MAT_lcov_unobs, iunobs, icol, 0.0); 	// ... set upper triangular matrix to zero 
		}
	}			
	
	//... mat_L is the lower triangular matrix with real and positive diagonal entries from Cholesky decomposition of Variance matrix: Sigma = L * L'
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
		clvar_unobs[iunobs] = gsl_matrix_get( MAT_lcov_unobs, iunobs, iunobs); 
		
	}
	icount =0; 
	for (int iunobs=1; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol =0; icol < iunobs ; icol++){
			clcov_unobs[icount] = gsl_matrix_get( MAT_lcov_unobs, iunobs, icol); 
			icount++; 
		}
	}
 
    
	for(int im = 0; im < NUM_measures_thetac; im++){

		cmu_measures_thetac[im] = temp_x[i];		i++; 
		
		calpha_measures_thetac[im] = temp_x[i];		i++; 

		csigma_error_thetac[im] = temp_x[i];		i++; 
	}
	
	for(int im = 0; im < NUM_measures_thetan; im++){

		cmu_measures_thetan[im]  = temp_x[i];		i++; 
		
		calpha_measures_thetan[im]  = temp_x[i];	i++; 

		csigma_error_thetan[im]  = temp_x[i];		i++; 
		
	}
	
	
	//... controls in the measurements
	for(int im = 0; im < NUM_measures_thetac; im++){
		// covariance of thetac 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			cbeta_measures_thetac[im][ic] = temp_x[i];		i++; 
		}
	}
	
	for(int im = 0; im < NUM_measures_thetan; im++){
		// covariance of thetac 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			cbeta_measures_thetan[im][ic] = temp_x[i];		i++; 
		}
	}

		
	for (int icat =0; icat < NUM_icat_unobs; icat++){
		
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
			cmu_unobs[icat][iunobs]  = temp_x[i];		i++; 
		}
	}
	
	
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){			
		cmu_unobs0[iunobs] = temp_x[i];		i++; 
	}

}


void printpar_loglike(const char *fname){

	FILE * outfile; outfile  = fopen (fname, "w"); 
	
    int i = 0;

	for(int im = 2; im < NUM_measures_thetan; im++){
 	fprintf(outfile, "%.8f \n",   cutoff_measure_thetan[im-2]);		i++; 
 	}
 		
	// unconditional distribution of unobservables
		for (int iunobs=0; iunobs < NUM_unobs_loglike; iunobs++){
			fprintf(outfile, "%.8f \n",   gsl_matrix_get(MAT_cov_unobs, iunobs, iunobs) );	i++; 
		}
		for (int iunobs=1; iunobs < NUM_unobs_loglike; iunobs++){
		for (int icol =0; icol < iunobs ; icol++){	
			fprintf(outfile, "%.8f \n",   gsl_matrix_get(MAT_cov_unobs, iunobs, icol) );	i++; 
		}}
		
    	
	for(int im = 0; im < NUM_measures_thetac; im++){

		fprintf(outfile, "%.8f \n",   cmu_measures_thetac[im] );		i++; 
		
		fprintf(outfile, "%.8f \n",   calpha_measures_thetac[im] );		i++; 

		fprintf(outfile, "%.8f \n",   csigma_error_thetac[im] );		i++; 
	}
	
	for(int im = 0; im < NUM_measures_thetan; im++){

		fprintf(outfile, "%.8f \n",   cmu_measures_thetan[im] );		i++; 
		
		fprintf(outfile, "%.8f \n",   calpha_measures_thetan[im] );		i++; 

		fprintf(outfile, "%.8f \n",   csigma_error_thetan[im] );		i++; 
	}

	
	//... controls in the measurements
	for(int im = 0; im < NUM_measures_thetac; im++){
		// covariance of thetac 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			fprintf(outfile, "%.8f \n",   cbeta_measures_thetac[im][ic] );	i++; 
		}
	}
	
	for(int im = 0; im < NUM_measures_thetan; im++){
		// covariance of thetac 
		for (int ic = 0; ic < NUM_measures_control; ic++){
			fprintf(outfile, "%.8f \n",   cbeta_measures_thetan[im][ic] );	i++; 
		}
	}

	
	for (int icat =0; icat < NUM_icat_unobs; icat++){
		
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
			fprintf(outfile, "%.8f \n",   cmu_unobs[icat][iunobs] );		i++; 
		}
	}
	
	
	
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){			
		fprintf(outfile, "%.8f \n",   cmu_unobs0[iunobs] );		i++; 
	}

	fclose(outfile); 
	
}


int rmvnorm(const gsl_rng *r, const int n, const gsl_vector *mean, const gsl_matrix *var, gsl_vector *result){
/* multivariate normal distribution random number generator */
/*
*	n	dimension of the random vetor
*	mean	vector of means of size n
*	var	variance matrix of dimension n x n
*	result	output variable with a sigle random vector normal distribution generation
*/
int k;
gsl_matrix *work = gsl_matrix_alloc(n,n);

gsl_matrix_memcpy(work,var);
gsl_linalg_cholesky_decomp(work);

for(k=0; k<n; k++)
	gsl_vector_set( result, k, gsl_ran_ugaussian(r) );

gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result,mean);

gsl_matrix_free(work);

return 0;
}


int rmvnorm_logpdf(const int M, const gsl_vector * x, const gsl_vector * mu, const gsl_matrix * L, double * result){
/* Compute the log of the probability density function at a given quantile
 * vector for a multivariate Gaussian distribution using the Cholesky
 * decomposition of the variance-covariance matrix.
 * x       vector of quantiles (dimension d)
 * mu      mean vector (dimension d)
 * L       matrix resulting from the Cholesky decomposition of
 *         variance-covariance matrix Sigma = L L^T (dimension d x d)
 * result  output of the density (dimension 1)
 * work    vector used for intermediate computations (dimension d)
 */
 
 
 	// work stores prediction error 
 	gsl_vector *work = gsl_vector_alloc(M); 
 
      double quadForm;        /* (x - mu)' Sigma^{-1} (x - mu) */
      double logSqrtDetSigma; /* log [ sqrt(|Sigma|) ] */

      /* compute: work = x - mu */
      for (int i = 0; i < M; ++i)
        {
          double xi = gsl_vector_get(x, i);
          double mui = gsl_vector_get(mu, i);
          gsl_vector_set(work, i, xi - mui);
        }

      /* compute: work = L^{-1} * (x - mu) */
      gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, L, work);

      /* compute: quadForm = (x - mu)' Sigma^{-1} (x - mu) */
      gsl_blas_ddot(work, work, &quadForm);

      /* compute: log [ sqrt(|Sigma|) ] = sum_i log L_{ii} */
      logSqrtDetSigma = 0.0;
      for (int i = 0; i < M; ++i)
        {
          double Lii = gsl_matrix_get(L, i, i);
          logSqrtDetSigma += log(Lii);
        }

  *result = -0.5*quadForm - logSqrtDetSigma - 0.5*M*log(2.0*M_PI);
  
  gsl_vector_free(work); 
      
return 0;  

}


// ================================================================


double density_measures_thetac(int ifirm, double thetac){
	
	double small_im = 1.0e-30; 

	double fc=1.0; 
	
 	 	//... measurement for thetac 
		for(int im = 0; im < NUM_measures_thetac; im++){

			double Xbeta = cmu_measures_thetac[im]; 
			
			for (int ix =0; ix < NUM_measures_control; ix++) Xbeta += DATA_init_measures_control[ix][ifirm] * cbeta_measures_thetac[im][ix]; 
			
			double xb = Xbeta + calpha_measures_thetac[im] * thetac;  
			
			double f_im = gsl_ran_gaussian_pdf( DATA_init_measures_thetac[im][ifirm] - xb ,  csigma_error_thetac[im] );
			if (f_im < small_im ) f_im = small_im; 
						
			fc *= f_im;   

		}
	
//	if (fc<1.0e-60) fc = 1.0e-60; 
	
	return fc; 
}


double density_measures_thetan(int ifirm, double thetan){
// Note: DATA_init_measures_thetan[im][ifirm] is a double variable 
	
	double small_im = 1.0e-30; 

	
	double fn=1.0; 
	
 	 	//... measurement for thetan 	
		for(int im = 0; im < 2; im++){
			
			double Xbeta = cmu_measures_thetan[im]; 
			
			for (int ix =0; ix < NUM_measures_control; ix++) Xbeta += DATA_init_measures_control[ix][ifirm] * cbeta_measures_thetan[im][ix]; 
			
			double xb = Xbeta + calpha_measures_thetan[im] * thetan;  
			
			// Prob(eps + xb > 0) = Prob( eps > -xb ) = Prob ( eps < xb ) 
			double pr1, pr0; 
	     	pr1 = gsl_cdf_gaussian_P(xb, csigma_error_thetan[im]); 
			pr0 = 1.0 - pr1; 
			
			double f_im = DATA_init_measures_thetan[im][ifirm] * pr1 + (1.0-DATA_init_measures_thetan[im][ifirm])*pr0;
			if (f_im < small_im ) f_im = small_im; 
						
			fn *= f_im; 
				
		} 	
	
		for(int im = 2; im < NUM_measures_thetan; im++){
			
			double Xbeta = cmu_measures_thetan[im]; 
			
			for (int ix =0; ix < NUM_measures_control; ix++) Xbeta += DATA_init_measures_control[ix][ifirm] * cbeta_measures_thetan[im][ix]; 
			
			double xb = Xbeta + calpha_measures_thetan[im] * thetan;  
			
			double pr[3]; 
			pr[0] = gsl_cdf_gaussian_P(0.0                         - xb, csigma_error_thetan[im]);                 // Prob ( eps + xb < cutoff0 )
			pr[1] = gsl_cdf_gaussian_P(cutoff_measure_thetan[im-2] - xb, csigma_error_thetan[im]) - pr[0];         // Prob ( cutoff0 < eps + xb < cutoff1 )
	     	pr[2] = 1.0 - pr[0] - pr[1]; 
	     	
	     	int ipr =  floor( DATA_init_measures_thetan[im][ifirm] );  // ranged: 0.0, 1.0, 2.0 
	     	
			double f_im = pr[ipr] ; 			
			
			if (f_im < small_im ) f_im = small_im; 
			if (f_im >= 1.0-small_im ) f_im = 1.0-small_im; 			
			
			fn *= f_im;   
	
		} 	

//	if (fn<1.0e-60) fn = 1.0e-60; 

	return fn; 	

}


double density(int ifirm){
	
	#ifdef DEBUG
		if (NUM_icat_unobs != 1) {
			cerr<< "\n Errors in logdensity: NUM_icat_unobs " << NUM_icat_unobs << " != 1 \n Exiting... \n "; 
			exit(1); 
		}
		if (NUM_measures_control != 3){
			cerr<< "\n Errors in logdensity: NUM_measures_control " << NUM_measures_control << " != 3 \n Exiting... \n "; 
			exit(1); 
		}
	#endif 
	
	
	double icat_control[NUM_icat_unobs]; 		
	for (int icat=0; icat < NUM_icat_unobs; icat++) icat_control[icat] = DATA_init_icat_control[icat][ifirm]; 
			
	double mu[NUM_unobs_loglike]; 
	for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
		mu[iunobs] = cmu_unobs0[iunobs]; 	
		for (int icat=0; icat < NUM_icat_unobs; icat++) mu[iunobs] += cmu_unobs[icat][iunobs]*icat_control[icat]; 
	}
	
	
	gsl_vector * vec_unobs  = gsl_vector_alloc(NUM_unobs_loglike); 	
	
		
	double ave_prob = 0.0; 
	for (int issp=0; issp < NUM_simsp_loglike; issp++){
 	 	
 	 	//... simulate multivariate normal distribution ... 	
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
 	 		gsl_vector_set( vec_unobs, iunobs, simp_rng_unobs_loglike[issp][iunobs]); 
 	 	}
 	 	gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, MAT_lcov_unobs, vec_unobs );
		
		//... add mean 
		double thetac = gsl_vector_get(vec_unobs, 0) + mu[0];
		double thetan = gsl_vector_get(vec_unobs, 1) + mu[1];
		
 	 	double fc = 1.0;  	 	double fn = 1.0;  	 	 	 	
 	 	
		fc = density_measures_thetac(ifirm, thetac); 
		 	
	 	fn = density_measures_thetan(ifirm, thetan); 
 	 	
		ave_prob += (fc*fn); 

 	}
		
	ave_prob /= NUM_simsp_loglike; 
	
	gsl_vector_free( vec_unobs ); 
		
	return ave_prob; 
	
}


double logdensity(int ifirm){

//	double smallpos = 1.0e-30; 	
		
	double den = density(ifirm); 
// 	if (den < smallpos) den = smallpos;
 	
	return log(den); 		
	
}


double loglikelihood(const gsl_vector *x, void *void_params){

	getpar_loglike(x, Nvar_loglike); // get the structure parameters from the input vector
	
	double sumlog = 0.0; 
	
	
	for (int ifirm=0; ifirm<NUM_loglike_ids; ifirm++){
		sumlog += logdensity(ifirm); 
	}
		
	return (-sumlog); 
	
}


double pred_measure(int num_control, double *measures_control, double theta, double calpha,double cmu, double *cbeta){
			
			double Xbeta = cmu; 
			
			for (int ix =0; ix < num_control; ix++) Xbeta += measures_control[ix] * cbeta[ix]; 
			
			double xb = Xbeta + calpha * theta;  
		
		return xb; 		
}


double pred_measure_dummy(int num_control, double *measures_control, double theta, double cmu, double *cbeta, double calpha, double csigma_error){
			
			double Xbeta = cmu; 
			
			for (int ix =0; ix < num_control; ix++) Xbeta += measures_control[ix] * cbeta[ix]; 
			
			double xb = Xbeta + calpha * theta;  
			
	     	double mean_pr1 = gsl_cdf_gaussian_P(xb, csigma_error); 
		
		return mean_pr1; 		
}



void test_loglike(){

	if (NUM_DATA_init != NUM_loglike_ids){
		printf("\n Errors in number of ids in loglike: NUM_DATA_init %4d != NUM_loglike_ids %4d \n ... exiting ... \n", NUM_DATA_init, NUM_loglike_ids); 
	}
		
	const int Nx = Nvar_loglike; 
	gsl_vector *x = gsl_vector_alloc(Nx); 
	
	getvec_loglike(x, Nvar_loglike); 
	getpar_loglike(x, Nvar_loglike); 
		
    double temp_dis = loglikelihood(x,x); 
		
    printf ("\n\n Loglikelihood function, value = %f \n", temp_dis);
		
	gsl_vector_free(x); 

}


double predict_factor(const double mu_prior, const double var_prior, const int num_measures, double y_measures[], double calpha_measures[], double csigma_error[], double cmu_measures[], double cbeta_measures[][NUM_measures_control], double measures_control[NUM_measures_control]){
/*
If the measurement system is linear, the posterior distribution is characterized by Kalman filter. 
See Handbook of Economic Forecasting, Chapter 7 by Harvey, and Cunha & Heckman (2006 JHR) 
*/
	
	gsl_vector * vec_invresid = gsl_vector_alloc(num_measures); 
	gsl_vector * vec_resid = gsl_vector_alloc(num_measures);  
	gsl_matrix * var_resid = gsl_matrix_alloc(num_measures, num_measures);
	gsl_matrix * var_error = gsl_matrix_alloc(num_measures, num_measures);

	// ... matrix for factor loadings, degenerated to a vector when there is only one factor ... 	
	gsl_vector * coef_alpha = gsl_vector_alloc(num_measures); 
		
	gsl_vector_set_zero(vec_resid); gsl_vector_set_zero(vec_invresid); 
	gsl_matrix_set_zero(var_resid); gsl_matrix_set_zero(var_error); 
	gsl_vector_set_zero(coef_alpha); 
	
		
 	 	//... measurement for thetac 
		for(int im = 0; im < num_measures; im++){
			
			double Xbeta = cmu_measures[im]; 
			
			for (int ix =0; ix < NUM_measures_control; ix++) Xbeta += measures_control[ix] * cbeta_measures[im][ix]; 
			
			double residual_y = y_measures[im] - Xbeta - calpha_measures[im] * mu_prior; 
			
			
			gsl_vector_set(vec_resid, im, residual_y); 			
			
			gsl_matrix_set(var_error, im, im, csigma_error[im]*csigma_error[im]); 
			
			gsl_vector_set(coef_alpha, im, calpha_measures[im]); 
						
		}
	
	// variance of forecasting error: var_resid = (coef_alpha * var_prior * coef_alpha' + var_error ),  var_prior is a number  
	gsl_matrix_set_zero(var_resid); 
	gsl_blas_dger (var_prior, coef_alpha, coef_alpha, var_resid); // These functions compute the rank-1 update A = \alpha x y^T + A of the matrix A.
	//gsl_blas_dgemm (CblasNoTrans, CblasTrans,   1.0, var_prior, coef_alpha, 0.0, var_resid);  // These functions compute sum C = \alpha op(A) op(B) + \beta C
	//gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, coef_alpha, var_resid, 0.0, var_resid); 
	gsl_matrix_add (var_resid, var_error);  // This function adds the elements of matrix b to the elements of matrix a
	
	
	// double thetac = mu_prior + (var_prior * coef_alpha') * (var_resid^(-1) * vec_resid); 
	
	// ... inverse of var_resid 
	gsl_linalg_cholesky_decomp(var_resid); 
	gsl_linalg_cholesky_invert(var_resid); 
	
	// ... vec_invresid = var_resid^(-1) * vec_resid
	gsl_blas_dgemv (CblasNoTrans, 1.0, var_resid, vec_resid, 0.0, vec_invresid); // compute sum y = \alpha op(A) x + \beta y
	
	
	// double thetac = mu_prior + (var_prior * coef_alpha') * vec_invresid, where vec_invresid = (var_resid^(-1) * vec_resid); 
	double theta = mu_prior; 
	for(int im = 0; im < num_measures; im++){ 
		
		theta += var_prior * gsl_vector_get(coef_alpha, im) * gsl_vector_get(vec_invresid, im); 
	}
	
	
	gsl_vector_free(vec_invresid);
	gsl_vector_free(vec_resid);
	gsl_matrix_free(var_resid);
	gsl_matrix_free(var_error);
	gsl_vector_free(coef_alpha);
	
	return theta; 
	
}


double predict_thetac_logdensity(double thetac, void *params){
			
	int ifirm = *(int *) params;
		
	int iunobs = 0; 
	double mu = cmu_unobs0[iunobs]; 	
	for (int icat=0; icat < NUM_icat_unobs; icat++) mu += cmu_unobs[icat][iunobs]*DATA_init_icat_control[icat][ifirm]; 
		
	double pfc = ( gsl_ran_gaussian_pdf( thetac - mu ,  1.0 ) ); 
		 	
	double fc = density_measures_thetac(ifirm, thetac); 

	if (pfc < 1.0e-30) pfc = 1.0e-30;   
	if (fc < 1.0e-30) fc = 1.0e-30;   
	 	 	
 	double f = pfc*fc; 
		
	return ( - log(f)); 

}


double predict_thetac(int ifirm){

		double x_r = 0.0; 
    	    	
		double x_lo = -10.0; 
		double x_hi =  10.0; 

		
		x_r = (x_lo + x_hi)/2.0; 
    	
				
		int status, iter, max_iter;
		 	   			
			gsl_function F; 
       		F.function = &predict_thetac_logdensity;
       		F.params   = &ifirm;
			       					
			const gsl_min_fminimizer_type *T;
  			gsl_min_fminimizer *s;
  			
		    T = gsl_min_fminimizer_brent;
  		    s = gsl_min_fminimizer_alloc (T);
  		    gsl_min_fminimizer_set (s, &F, x_r, x_lo, x_hi);
			
			iter = 0; 
			max_iter = 50; 
			
			double epsrel = 0.0; 			
			double epsabs = 1.0e-3;
			
  			do
    		{
      			iter++;
      			status = gsl_min_fminimizer_iterate (s);

      			x_r  = gsl_min_fminimizer_x_minimum (s);
      			x_lo = gsl_min_fminimizer_x_lower (s);
      			x_hi = gsl_min_fminimizer_x_upper (s);

				status = gsl_min_test_interval (x_lo, x_hi, epsabs, epsrel);
	  
    		}
  		while (status == GSL_CONTINUE && iter < max_iter);
	
	gsl_min_fminimizer_free (s);
	 
	double thetac = x_r; 
	
	return thetac; 
	
}


double predict_thetan_logdensity(double thetan, void *params){
			
	int ifirm = *(int *) params;
		
	int iunobs = 1; 
	double mu = cmu_unobs0[iunobs]; 	
	for (int icat=0; icat < NUM_icat_unobs; icat++) mu += cmu_unobs[icat][iunobs]*DATA_init_icat_control[icat][ifirm]; 
		
	double pfn = gsl_ran_gaussian_pdf( thetan - mu ,  1.0 ); 
		
	double fn = density_measures_thetan(ifirm, thetan); 
	 	
	if (pfn < 1.0e-60) pfn = 1.0e-60;   
	if (fn < 1.0e-60) fn = 1.0e-60;   
		
	double f = pfn * fn; 
	 		
	return ( - log(f)); 

}


double predict_thetan(const int ifirm0){
		
		int ifirm = ifirm0; 
		
		double x_r = 0.0; 
    	    	
		double x_lo = -10.0; 
		double x_hi =  10.0; 
		
		x_r = (x_lo + x_hi)/2.0; 
			
		int status, iter, max_iter;
		 	   			
			gsl_function F; 
       		F.function = &predict_thetan_logdensity;
       		F.params   = &ifirm;
			       					
			const gsl_min_fminimizer_type *T;
  			gsl_min_fminimizer *s;
  			
		    T = gsl_min_fminimizer_brent;
  		    s = gsl_min_fminimizer_alloc (T);
  		    gsl_min_fminimizer_set (s, &F, x_r, x_lo, x_hi);
			
			iter = 0; 
			max_iter = 50; 
			
			double epsrel = 0.0; 			
			double epsabs = 1.0e-3;
			
  			do
    		{
      			iter++;
      			status = gsl_min_fminimizer_iterate (s);

      			x_r  = gsl_min_fminimizer_x_minimum (s);
      			x_lo = gsl_min_fminimizer_x_lower (s);
      			x_hi = gsl_min_fminimizer_x_upper (s);

				status = gsl_min_test_interval (x_lo, x_hi, epsabs, epsrel);
			  
    		}
  		while (status == GSL_CONTINUE && iter < max_iter);
	
	gsl_min_fminimizer_free (s);
	 
	double thetan = x_r; 
	
	return thetan ; 	
}



void test_predict(){
// ... This program double check the predicted factors (thetac in particular). 
		
    int n; 
    ifstream infile;   
    
    n = NUM_DATA_init; 	   
	
	double thetac[n];  
	
	double sumsq_diff = 0.0; 
	double sumsq_resid = 0.0; 	
	double sumsq_resid_linear = 0.0; 	
				
	for (int ifirm=0; ifirm<NUM_loglike_ids ; ifirm++){  

		double mu_prior, var_prior; 
		
		double icat_control[NUM_icat_unobs]; 		
		for (int icat=0; icat < NUM_icat_unobs; icat++) icat_control[icat] = DATA_init_icat_control[icat][ifirm]; 
	
		double mu[NUM_unobs_loglike]; 
		for (int iunobs=0; iunobs< NUM_unobs_loglike; iunobs++){
			mu[iunobs] = cmu_unobs0[iunobs]; 	
			for (int icat=0; icat < NUM_icat_unobs; icat++) mu[iunobs] += cmu_unobs[icat][iunobs]*icat_control[icat]; 
		}
		
				
		double measures_control[NUM_measures_control];		
		for (int ix =0; ix < NUM_measures_control; ix++) {
			measures_control[ix]       = DATA_init_measures_control[ix][ifirm]; 
		}
		
		
		// ... predict thetac 
		mu_prior = mu[0]; 	var_prior = 1.0; 
					
		thetac[ifirm] = predict_thetac(ifirm); 
		double wgt    = gsl_ran_gaussian_pdf( thetac[ifirm] - mu_prior, var_prior);

		
		// ... predict thetac based on Kalman filter, yielding to the same results as maximizing the loglikelhood because the measures are linear 
		double thetac_measures[NUM_measures_thetac];		
		for(int im = 0; im < NUM_measures_thetac; im++){
			thetac_measures[im] = DATA_init_measures_thetac[im][ifirm]; 
		}	
						
		double thetac_linear = predict_factor(mu_prior, var_prior, NUM_measures_thetac, thetac_measures, calpha_measures_thetac, csigma_error_thetac, cmu_measures_thetac, cbeta_measures_thetac, measures_control); 
		double wgt_linear    = gsl_ran_gaussian_pdf( thetac_linear - mu_prior, var_prior);

		
		// ... check predicted thetac is the same as the the prediction based on Kalman filter since the measurements are linear
		sumsq_diff += ( (thetac[ifirm] - thetac_linear) * (thetac[ifirm] - thetac_linear) ); 

		
		// check if predicted thetac minimized the sum of forecasting error in measurements weighted by its prior pdf 
		for(int im = 0; im < NUM_measures_thetac; im++){

			double Xbeta = cmu_measures_thetac[im]; 
			for (int ix =0; ix < NUM_measures_control; ix++) Xbeta += DATA_init_measures_control[ix][ifirm] * cbeta_measures_thetac[im][ix]; 
				
			double residual_y = DATA_init_measures_thetac[im][ifirm] - pred_measure(NUM_measures_control, measures_control, thetac[ifirm], calpha_measures_thetac[im], cmu_measures_thetac[im], cbeta_measures_thetac[im]); 
			
			sumsq_resid += residual_y * residual_y * wgt; 	
			
			
			double residual_y_linear = DATA_init_measures_thetac[im][ifirm] - Xbeta - calpha_measures_thetac[im] * thetac_linear; 
			
			sumsq_resid_linear += residual_y_linear * residual_y_linear * wgt_linear; 	
				
		}
			
	}
	
	
	printf("\n\t test_predict(): sumsq_diff  %f, \t sumsq_resid %f : sumsq_resid_linear %f \n\n", sumsq_diff, sumsq_resid, sumsq_resid_linear); 
	
	if ( sumsq_diff  > 1.0e-4 || fabs (sumsq_resid - sumsq_resid_linear) > 1.0e-4) {printf("\n Error!!! Please recheck the predicted thetas \n exiting ...\n"); exit(1); }

}

void predict_unobs(){

    int n; 
    ifstream infile;   	
	
	FILE * out_smmfile; out_smmfile  = fopen ("out_data_new/predicted_thetas.txt", "w"); 
	
	n = NUM_DATA_init; 	   	
	
	double thetac[n], thetan[n]; 
	
	for (int ids=0; ids<NUM_loglike_ids ; ids++){ 
				
		thetac[ids] = predict_thetac(ids);
		
		thetan[ids] = predict_thetan(ids);		
		
		int idnew = ids+1; 
	    fprintf(out_smmfile, "%d \t %f \t %f \n", idnew, thetac[ids], thetan[ids]) ; 
							
	}
	
	fclose(out_smmfile); 
	
	printf("\n\t predict_unobs(): printing out_data_new/predicted_thetas.txt \n "); 
	
}


void estiamte_loglike(){

	ifstream infile; 	
	infile.open("estpar_loglike_iter.txt"); readpar_loglike(500, infile);	infile.close(); 	
	
	
	const int Nx = Nvar_loglike; 
	gsl_vector *x = gsl_vector_alloc(Nx); 

    double init_size      = Est_stepsize_loglike; 
    double size_tolerance = Est_tol_loglike; 
    double max_iter       = Est_maxiter_loglike; 
	
	
		
	getvec_loglike(x, Nvar_loglike); 
	getpar_loglike(x, Nvar_loglike); 

	
	for (int i_s =0; i_s<Nx; i_s++) printf(" \t x [%3d] = %5.6f ; \n", i_s, gsl_vector_get(x,i_s));
	
		
    // !!!!!!!!!!! PARAMETER ESTIMATION STARTS !!!!!!!!!!	
    time_t start0,end0;	double timediff; 
	
    time (&start0);
	
    double temp_dis = loglikelihood(x,x); 
		
    time (&end0);	timediff = difftime (end0,start0);
	
    printf ("\n\n It takes %f minutes to evaluate loglikelihood function, Initial value = %f \n", timediff/60, temp_dis);
		
	
    // Start estimation: 
    printf ("\n\n Estimation Starts!!! \n");
	
    time_t start_t, end_t;	double diff_t; 
    time (&start_t);
    
    gsl_vector *opt_x = gsl_vector_alloc(Nx);	
    gsl_vector *opt_f = gsl_vector_alloc(2); 
    
    
    {
	
    gsl_vector * x_iter = gsl_vector_alloc(Nx);
	
    const gsl_multimin_fminimizer_type *T =  gsl_multimin_fminimizer_nmsimplex;   // gsl_multimin_fminimizer_nmsimplex2;
    gsl_multimin_fminimizer *s = NULL;
    gsl_multimin_function minex_func;
   	
    // Initialize method and iterate
    minex_func.f = &loglikelihood;
    minex_func.n = Nx;
    minex_func.params = x;
	
		
    // Set initial step sizes according to init_size   
    gsl_vector *ss = gsl_vector_alloc (Nx);
    gsl_vector_set_all (ss, init_size);

    s = gsl_multimin_fminimizer_alloc (T, Nx);
    gsl_multimin_fminimizer_set (s, &minex_func, x, ss);

    int iter = -1;    int status;    double size;		
    do {

	iter++;
	status = gsl_multimin_fminimizer_iterate(s);
	
	if (status) break;
   
	size = gsl_multimin_fminimizer_size (s);
	status = gsl_multimin_test_size (size, size_tolerance);
	
			         
	if ( iter%200==0 || status == GSL_SUCCESS){
		
		#ifdef DEBUG 
		printf ("%5d f() = %18.4f,  size = %1.5f \n", iter, s->fval, size);
		for (int i_s =0; i_s<Nx; i_s++)
	    	printf(" \t x [%3d] = %5.6f ; \n", i_s, gsl_vector_get (s->x, i_s));
		#endif
		
			
	    for (int i_s =0; i_s<Nx; i_s++){
		gsl_vector_set(x_iter, i_s, gsl_vector_get (s->x, i_s)); 
	    }
			
	    getpar_loglike(x_iter, Nvar_loglike); 
	    printpar_loglike("estpar_loglike_iter.txt");	
	
	}
	
    } while (status == GSL_CONTINUE && iter < max_iter);

    
    // store the optimal points;
    for (int i=0; i<Nx; i++)
	gsl_vector_set(opt_x, i, gsl_vector_get (s->x, i)); 	
	
    // store the function value at optimal point and the step size
    gsl_vector_set(opt_f, 0, (s->fval)); 
    gsl_vector_set(opt_f, 1, gsl_multimin_fminimizer_size (s) ); 
	
    gsl_vector_free(ss);  gsl_multimin_fminimizer_free (s); 
    gsl_vector_free(x_iter); 
	
    }
    
    
    time (&end_t);	diff_t = difftime (end_t, start_t);
    printf ("\n It took you %f hours to finish Estimation: f %f, size %f \n", diff_t/3600, gsl_vector_get(opt_f,0), gsl_vector_get(opt_f, 1) );
		
	test_loglike();  predict_unobs(); 	test_predict(); 
	
		
	gsl_vector_free(x); 
			
}

